
rm(list=ls())
library(metafor)
###################
#### Loading Data and AK Estimators
###################
PATH="C:/Users/hong_/Desktop/World Development/FinalCodes_Data/"
DT<-as.data.frame(read.csv(paste(PATH, 'data_sorted.csv', sep='')))
source(paste(PATH, 'AK_Estimators.R', sep=''))


###################
#### Creating Dataset and attache to R working directory
###################
RegData <- as.data.frame(cbind(effect=DT$pcc, se=DT$sepcc, DT))
attach(RegData)


Table06 <- as.data.frame(matrix(NA, nrow=48, ncol=8))
colnames(Table06) <- c('(1) FE', '(2) FEwSE', '(3) RE', '(4) REwSE', '(5) AK1', '(6) AK1wSE', '(7) AK2', '(8) AK2wSE')
rownames(Table06) <- c('intrcpt', 'se(intrcpt)', 'sepcc', 'se(sepcc)', 'fdilog', 'se(fdilog)', 'fdiother',
                       'se(fdiother)', 'entretype2', 'se(entretype2)', 'entretype3', 'se(entretype3)', 'entreother',
                       'se(entreother)', 'oecd', 'se(oecd)',
                       'nonoecd', 'se(nonoecd)', 'crosssection', 'se(crosssection)', 'panelfe', 'se(panelfe)',
                       'lagdv', 'se(lagdv)', 'busenviron', 'se(busenviron)', 'iv', 'se(iv)','journal',  
                       'se(journal)', '__', '___', 'tau', 'Beta1', 'se(Beta1)', 'Beta2', 'se(Beta2)', 'Beta3', 'se(Beta3)', 'obs', 'parameter',
                       'LLK', 'AIC', 'BIC', 'Test1', 'p-valueT1', 'Test2', 'p-valueT2')


###################
#### Making Table 6 - Upper Panel (Estimation Results)
###################
#Regression 01 - FE
regFE <- rma.uni(pcc ~ 1 + fdilog + fdiother + entretype2 + entretype3 + entreother + 
                   oecd + nonoecd + crosssection + panelfe + lagdv + busenviron + iv
                 + journal, vi=sepcc*sepcc, method='FE')
regFE_Clustered <- robust(regFE, id)
b <- regFE_Clustered$b
regFE_Clustered <- cbind(regFE_Clustered$b, regFE_Clustered$se)
FE_Estimates_Vector <- as.data.frame(matrix(0, nrow=(nrow(regFE_Clustered)*2), ncol=1))
FE_Estimates_Vector[seq(1,nrow(FE_Estimates_Vector),2),1] <- regFE_Clustered[,1]
FE_Estimates_Vector[seq(2,nrow(FE_Estimates_Vector),2),1] <- regFE_Clustered[,2]
FE_Estimates_Vector<-c(FE_Estimates_Vector[c(1:2),1], NA, NA, FE_Estimates_Vector[c(3:nrow(FE_Estimates_Vector)),1])
regFE<-summary(regFE)
Table06[c(1:30),1] <- FE_Estimates_Vector
Table06[c(40:44),1] <-  c(length(pcc), regFE$parms, regFE$fit.stats[1,1], regFE$fit.stats[3,1], regFE$fit.stats[4,1])
Pvalue_regFE <- as.matrix(pt(regFE_Clustered[,1]/regFE_Clustered[,2], df=(length(pcc)-regFE$parms)))
Pvalue_regFE <- cbind(coef=b, pval=Pvalue_regFE, sig=(1*(Pvalue_regFE<=0.01) + 1*(Pvalue_regFE<=0.05) + 1*(Pvalue_regFE<=0.10)))



#Regression 02 - FEwSE
regFEwSE <- rma.uni(pcc ~ 1 + sepcc + fdilog + fdiother + entretype2 + entretype3 + entreother 
                    + oecd + nonoecd + crosssection + panelfe + lagdv + busenviron + iv
                    + journal, vi=sepcc*sepcc, method='FE')
regFEwSE_Clustered <- robust(regFEwSE, id)
b <- regFEwSE_Clustered$b
regFEwSE_Clustered <- cbind(regFEwSE_Clustered$b, regFEwSE_Clustered$se)
FEwSE_Estimates_Vector <- as.data.frame(matrix(0, nrow=(nrow(regFEwSE_Clustered)*2), ncol=1))
FEwSE_Estimates_Vector[seq(1,nrow(FEwSE_Estimates_Vector),2),1] <- regFEwSE_Clustered[,1]
FEwSE_Estimates_Vector[seq(2,nrow(FEwSE_Estimates_Vector),2),1] <- regFEwSE_Clustered[,2]
regFEwSE<-summary(regFEwSE)
Table06[c(1:30),2] <- FEwSE_Estimates_Vector
Table06[c(40:44),2] <-  c(length(pcc), regFEwSE$parms, regFEwSE$fit.stats[1,1], regFEwSE$fit.stats[3,1], regFEwSE$fit.stats[4,1])
Pvalue_regFEwSE <- as.matrix(pt(regFEwSE_Clustered[,1]/regFEwSE_Clustered[,2], df=(length(pcc)-regFEwSE$parms)))
Pvalue_regFEwSE <- cbind(coef=b, pval=Pvalue_regFEwSE, sig=(1*(Pvalue_regFEwSE<=0.01) + 1*(Pvalue_regFEwSE<=0.05) + 1*(Pvalue_regFEwSE<=0.10)))




#Regression 03 - RE
regRE <- rma.uni(pcc ~ 1 + fdilog + fdiother + entretype2 + entretype3 + entreother 
                 + oecd + nonoecd + crosssection + panelfe + lagdv + busenviron + iv
                 + journal, vi=sepcc*sepcc, method='ML')
regRE_Clustered <- robust(regRE, id)
b <- regRE_Clustered$b
regRE_Clustered <- cbind(regRE_Clustered$b, regRE_Clustered$se)
RE_Estimates_Vector <- as.data.frame(matrix(0, nrow=(nrow(regRE_Clustered)*2), ncol=1))
RE_Estimates_Vector[seq(1,nrow(RE_Estimates_Vector),2),1] <- regRE_Clustered[,1]
RE_Estimates_Vector[seq(2,nrow(RE_Estimates_Vector),2),1] <- regRE_Clustered[,2]
RE_Estimates_Vector<-c(RE_Estimates_Vector[c(1:2),1], NA, NA, RE_Estimates_Vector[c(3:nrow(RE_Estimates_Vector)),1])
regRE<-summary(regRE)
Table06[c(1:30),3] <- RE_Estimates_Vector
Table06[c(33, 40:44),3]  <-  c(sqrt(regRE$tau2), length(pcc), regRE$parms, regRE$fit.stats[1,1], regRE$fit.stats[3,1], regRE$fit.stats[4,1])
Pvalue_regRE <- as.matrix(pt(regRE_Clustered[,1]/regRE_Clustered[,2], df=(length(pcc)-regRE$parms)))
Pvalue_regRE <- cbind(coef=b, pval=Pvalue_regRE, sig=(1*(Pvalue_regRE<=0.01) + 1*(Pvalue_regRE<=0.05) + 1*(Pvalue_regRE<=0.10)))




#Regression 04 - REwSE
regREwSE <- rma.uni(pcc ~ 1 + sepcc + fdilog + fdiother + entretype2 + entretype3 + entreother 
                    + oecd + nonoecd + crosssection + panelfe + lagdv + busenviron + iv
                    + journal, vi=sepcc*sepcc, method='ML')
regREwSE_Clustered <- robust(regREwSE, id)
b <- regREwSE_Clustered$b
regREwSE_Clustered <- cbind(regREwSE_Clustered$b, regREwSE_Clustered$se)
REwSE_Estimates_Vector <- as.data.frame(matrix(0, nrow=(nrow(regREwSE_Clustered)*2), ncol=1))
REwSE_Estimates_Vector[seq(1,nrow(REwSE_Estimates_Vector),2),1] <- regREwSE_Clustered[,1]
REwSE_Estimates_Vector[seq(2,nrow(REwSE_Estimates_Vector),2),1] <- regREwSE_Clustered[,2]
regREwSE<-summary(regREwSE)
Table06[c(1:30),4] <- REwSE_Estimates_Vector
Table06[c(33, 40:44),4] <-  c(sqrt(regREwSE$tau2), length(pcc), regREwSE$parms, regREwSE$fit.stats[1,1], regREwSE$fit.stats[3,1], regREwSE$fit.stats[4,1])
Pvalue_regREwSE <- as.matrix(pt(regREwSE_Clustered[,1]/regREwSE_Clustered[,2], df=(length(pcc)-regREwSE$parms)))
Pvalue_regREwSE <- cbind(coef=b, pval=Pvalue_regREwSE, sig=(1*(Pvalue_regREwSE<=0.01) + 1*(Pvalue_regREwSE<=0.05) + 1*(Pvalue_regREwSE<=0.10)))



#Regression 05 - AK1
AKdata <- as.data.frame(cbind(id, effect=pcc, se=sepcc, constant=1, 
                              fdilog, fdiother, entretype2, entretype3, entreother, 
                              oecd, nonoecd, crosssection, panelfe, lagdv,
                              busenviron, iv, journal))
AK1 <- AK1Estimator(AKdata)
regAK1_Clustered <- t(AK1$EstResults)[c(3:ncol(AK1$EstResults)), c(1,2)]
AK1_Estimates_Vector <- as.data.frame(matrix(0, nrow=(nrow(regAK1_Clustered)*2), ncol=1))
AK1_Estimates_Vector[seq(1,nrow(AK1_Estimates_Vector),2),1] <- regAK1_Clustered[,1]
AK1_Estimates_Vector[seq(2,nrow(AK1_Estimates_Vector),2),1] <- regAK1_Clustered[,2]
AK1_Estimates_Vector<-c(AK1_Estimates_Vector[c(1:2),1], NA, NA, AK1_Estimates_Vector[c(3:nrow(AK1_Estimates_Vector)),1])
Table06[c(1:30),5] <- AK1_Estimates_Vector
Table06[c(33:35, 40:44),5]  <-  c(AK1$EstResults[1,1], AK1$EstResults[1,2], AK1$EstResults[2,2], length(pcc), ncol(AK1$EstResults), AK1$LogLiklihood, AK1$aic, AK1$bic)
Pvalue_AK1 <- t(as.matrix(AK1$EstResults[c(1,3), 3:ncol(AK1$EstResults)]))
Pvalue_AK1 <- as.data.frame(cbind(Pvalue_AK1, sig=(1*(Pvalue_AK1[,2]<=0.01) + 1*(Pvalue_AK1[,2]<=0.05) + 1*(Pvalue_AK1[,2]<=0.10))))



#Regression 06 - AK1wSE
AKdata <- as.data.frame(cbind(id, effect=pcc, se=sepcc, constant=1, sepcc,
                              fdilog, fdiother, entretype2, entretype3, entreother, 
                              oecd, nonoecd, crosssection, panelfe, lagdv,
                              busenviron, iv, journal))
AK1wSE <- AK1Estimator(AKdata)
regAK1wSE_Clustered <- t(AK1wSE$EstResults)[c(3:ncol(AK1wSE$EstResults)), c(1,2)]
AK1wSE_Estimates_Vector <- as.data.frame(matrix(0, nrow=(nrow(regAK1wSE_Clustered)*2), ncol=1))
AK1wSE_Estimates_Vector[seq(1,nrow(AK1wSE_Estimates_Vector),2),1] <- regAK1wSE_Clustered[,1]
AK1wSE_Estimates_Vector[seq(2,nrow(AK1wSE_Estimates_Vector),2),1] <- regAK1wSE_Clustered[,2]
Table06[c(1:30),6] <- AK1wSE_Estimates_Vector
Table06[c(33:35, 40:44),6]  <-  c(AK1wSE$EstResults[1,1], AK1wSE$EstResults[1,2], AK1wSE$EstResults[2,2], length(pcc), ncol(AK1wSE$EstResults), AK1wSE$LogLiklihood, AK1wSE$aic, AK1wSE$bic)
Pvalue_AK1wSE <- t(as.matrix(AK1wSE$EstResults[c(1,3), 3:ncol(AK1wSE$EstResults)]))
Pvalue_AK1wSE <- cbind(Pvalue_AK1wSE, sig=(1*(Pvalue_AK1wSE[,2]<=0.01) + 1*(Pvalue_AK1wSE[,2]<=0.05) + 1*(Pvalue_AK1wSE[,2]<=0.10)))





#Regression 07 - AK2
AKdata <- as.data.frame(cbind(id, effect=pcc, se=sepcc, constant=1, 
                              fdilog, fdiother, entretype2, entretype3, entreother, 
                              oecd, nonoecd, crosssection, panelfe, lagdv,
                              busenviron, iv, journal))
AK2 <- AK2Estimator(AKdata)
regAK2_Clustered <- t(AK2$EstResults)[c(5:ncol(AK2$EstResults)), c(1,2)]
AK2_Estimates_Vector <- as.data.frame(matrix(0, nrow=(nrow(regAK2_Clustered)*2), ncol=1))
AK2_Estimates_Vector[seq(1,nrow(AK2_Estimates_Vector),2),1] <- regAK2_Clustered[,1]
AK2_Estimates_Vector[seq(2,nrow(AK2_Estimates_Vector),2),1] <- regAK2_Clustered[,2]
AK2_Estimates_Vector<-c(AK2_Estimates_Vector[c(1:2),1], NA, NA, AK2_Estimates_Vector[c(3:nrow(AK2_Estimates_Vector)),1])
Table06[c(1:30),7] <- AK2_Estimates_Vector
Table06[c(33:44),7]  <-  c(AK2$EstResults[1,1], 
                           AK2$EstResults[1,2], AK2$EstResults[2,2], 
                           AK2$EstResults[1,3], AK2$EstResults[2,3], 
                           AK2$EstResults[1,4], AK2$EstResults[2,4],
                           length(pcc), ncol(AK2$EstResults), AK2$LogLiklihood, AK2$aic, AK2$bic)
Pvalue_AK2 <- t(as.matrix(AK2$EstResults[c(1,3), 3:ncol(AK2$EstResults)]))
Pvalue_AK2 <- cbind(Pvalue_AK2, sig=(1*(Pvalue_AK2[,2]<=0.01) + 1*(Pvalue_AK2[,2]<=0.05) + 1*(Pvalue_AK2[,2]<=0.10)))





#Regression 08 - AK2wSE
AKdata <- as.data.frame(cbind(id, effect=pcc, se=sepcc, constant=1, sepcc,
                              fdilog, fdiother, entretype2, entretype3, entreother, 
                              oecd, nonoecd, crosssection, panelfe, lagdv,
                              busenviron, iv, journal))
AK2wSE <- AK2Estimator(AKdata)
regAK2wSE_Clustered <- t(AK2wSE$EstResults)[c(5:ncol(AK2wSE$EstResults)), c(1,2)]
AK2wSE_Estimates_Vector <- as.data.frame(matrix(0, nrow=(nrow(regAK2wSE_Clustered)*2), ncol=1))
AK2wSE_Estimates_Vector[seq(1,nrow(AK2wSE_Estimates_Vector),2),1] <- regAK2wSE_Clustered[,1]
AK2wSE_Estimates_Vector[seq(2,nrow(AK2wSE_Estimates_Vector),2),1] <- regAK2wSE_Clustered[,2]
Table06[c(1:30),8] <- AK2wSE_Estimates_Vector
Table06[c(33:44),8]  <-  c(AK2wSE$EstResults[1,1], 
                           AK2wSE$EstResults[1,2], AK2wSE$EstResults[2,2], 
                           AK2wSE$EstResults[1,3], AK2wSE$EstResults[2,3], 
                           AK2wSE$EstResults[1,4], AK2wSE$EstResults[2,4],
                           length(pcc), ncol(AK2wSE$EstResults), AK2wSE$LogLiklihood, AK2wSE$aic, AK2wSE$bic)
Pvalue_AK2wSE <- t(as.matrix(AK2wSE$EstResults[c(1,3), 3:ncol(AK2wSE$EstResults)]))
Pvalue_AK2wSE <- cbind(Pvalue_AK2wSE, sig=(1*(Pvalue_AK2wSE[,2]<=0.01) + 1*(Pvalue_AK2wSE[,2]<=0.05) + 1*(Pvalue_AK2wSE[,2]<=0.10)))








###################
#### Making Table 6 - Lower Panel (Test Results)
###################
para <- as.numeric(Table06[41,])
llik <- as.numeric(Table06[42,])
## Column 1 (FE)
TestPair1 <- c(2,1) # Column 2 (unrestricted) vs Column 1 (Restricted)
test1.p.value <- pchisq(2*(llik[TestPair1[1]]-llik[TestPair1[2]]), (para[TestPair1[1]]-para[TestPair1[2]]), lower.tail=FALSE)
TestPair2 <- c(3,1)
test2.p.value <- pchisq(2*(llik[TestPair2[1]]-llik[TestPair2[2]]), (para[TestPair2[1]]-para[TestPair2[2]]), lower.tail=FALSE)
Table06[45:48,1] <- c(TestPair1[1], test1.p.value, TestPair2[1], test2.p.value)

## Column 2 (FEwSE)
TestPair1 <- c(4,2) 
test1.p.value <- pchisq(2*(llik[TestPair1[1]]-llik[TestPair1[2]]), (para[TestPair1[1]]-para[TestPair1[2]]), lower.tail=FALSE)
Table06[45:46,2] <- c(TestPair1[1], test1.p.value)

## Column 3 (RE)
TestPair1 <- c(4,3)
test1.p.value <- pchisq(2*(llik[TestPair1[1]]-llik[TestPair1[2]]), (para[TestPair1[1]]-para[TestPair1[2]]), lower.tail=FALSE)
TestPair2 <- c(5,3)
test2.p.value <- pchisq(2*(llik[TestPair2[1]]-llik[TestPair2[2]]), (para[TestPair2[1]]-para[TestPair2[2]]), lower.tail=FALSE)
Table06[45:48,3] <- c(TestPair1[1], test1.p.value, TestPair2[1], test2.p.value)

## Column 4 (REwSE)
TestPair1 <- c(6,4)
test1.p.value <- pchisq(2*(llik[TestPair1[1]]-llik[TestPair1[2]]), (para[TestPair1[1]]-para[TestPair1[2]]), lower.tail=FALSE)
Table06[45:46,4] <- c(TestPair1[1], test1.p.value)

## Column 5 (AK1)
TestPair1 <- c(6,5)
test1.p.value <- pchisq(2*(llik[TestPair1[1]]-llik[TestPair1[2]]), (para[TestPair1[1]]-para[TestPair1[2]]), lower.tail=FALSE)
TestPair2 <- c(7,5)
test2.p.value <- pchisq(2*(llik[TestPair2[1]]-llik[TestPair2[2]]), (para[TestPair2[1]]-para[TestPair2[2]]), lower.tail=FALSE)
Table06[45:48,5] <- c(TestPair1[1], test1.p.value, TestPair2[1], test2.p.value)

## Column 6 (AK1wSE)
TestPair1 <- c(8,6)
test1.p.value <- pchisq(2*(llik[TestPair1[1]]-llik[TestPair1[2]]), (para[TestPair1[1]]-para[TestPair1[2]]), lower.tail=FALSE)
Table06[45:46,6] <- c(TestPair1[1], test1.p.value)

## Column 7 (AK2)
TestPair1 <- c(8,7)
test1.p.value <- pchisq(2*(llik[TestPair1[1]]-llik[TestPair1[2]]), (para[TestPair1[1]]-para[TestPair1[2]]), lower.tail=FALSE)
Table06[45:46,7] <- c(TestPair1[1], test1.p.value)


Table06 <- as.data.frame(Table06)
Table06[45,] <- colnames(Table06)[as.numeric( Table06[45,] )]
Table06[47,] <- colnames(Table06)[as.numeric( Table06[47,] )]

cat('Convergence Summary:\n',
    "AK1: ",  AK1$nlmOutput$message, '; # of iterations: ', AK1$nlmOutput$iterations,'\n',
    "AK1wSE: ",  AK1wSE$nlmOutput$message, '; # of iterations: ', AK1wSE$nlmOutput$iterations,'\n',
    "AK2: ",  AK2$nlmOutput$message, '; # of iterations: ', AK2$nlmOutput$iterations,'\n',
    "AK2wSE: ",  AK2wSE$nlmOutput$message, '; # of iterations: ', AK2wSE$nlmOutput$iterations,'\n')
detach(RegData)
View(Table06)
write.csv(as.data.frame(Table06), "EstimationResults/Table07.csv", row.names=TRUE)